Lagrangian method for multiple correlations in passive scalar advection 
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A Lagrangian method is introduced for calculating simultaneous n-point correlations of a passive 
scalar advected by a random velocity field, with random forcing and finite molecular diffusivity k. 
The method, which is here presented in detail, is particularly well suited for studying the k — > 
limit when the velocity field is not smooth. Efficient Monte Carlo simulations based on this method 
are applied to the Kraichnan model of passive scalar and lead to accurate determinations of the 
anomalous intermittency corrections in the fourth-order structure function as a function of the 
scaling exponent £ of the velocity field in two and three dimensions. Anomalous corrections are 
found to vanish in the limits £ — > and £ — > 2, as predicted by perturbation theory. 
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I. INTRODUCTION 



II. THE LAGRANGIAN METHOD 
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Robert Kraichnan's model of passive scalar advection 
by a white-in-time velocity field has been particularly fer- 
tile ground for theoreticians trying to develop a theory 
of intermittency jl)-^ (see also Refs. QH). Although the 
model leads to closed equations for multiple-point mo- 
ments, only second-order moments can be obtained in 
closed analytic form ||. Theoretical predictions differed 
as to the behavior of higher-order quantities, regarding in 
particular the survival or the vanishing of intermittency 
corrections (anomalies) in certain limits. Obtaining reli- 
able numerical results was thus an important challenge. 
Until recently numerical simulations have been based on 
the direct integration of the passive scalar partial differ- 
ential equation and have been limited to two dimensions 
Hfr|j|]. Such calculations are delicate; to wit, the dif- 
ficulty of observing for the second-order structure func- 
tion the known high-Pcclet number asymptotic scaling 
||. Also, the numerical scheme used in Refs. (5[0] in- 
volves a slightly anisotropic velocity field which is not 
expected to give exactly the right scaling laws for the 
passive scalar 0. 

Lagrangian methods for tackling the Kraichnan model 
and which require only the integration of ordinary dif- 
ferential equations were recently proposed independently 
by Frisch, Mazzino and Vergassola fnj and by Gat, Pro- 
caccia and Zeitak |tl| . Our goal here is to give a detailed 
presentation of the Lagrangian method and to present 
new results. 

In Section || we give the theoretical background of 
the Lagrangian method for general random velocity fields 
which need not be white-in-time. In Section III 



we inves- 



tigate the limit of vanishing molecular diffusivity which 
depends crucially on how nearby Lagrangian trajectories 
separate. We then turn to the Kraichnan model which 
is given a Lagrangian formulation (Section |Tv| ). Then, 
we show how it can be solved numerically by a Monte- 
Carlo method (Section 0) and present results in both 
two and three space dimensio ns (S ection VI). We make 
concluding remarks in Section VII . 



The (Eulerian) dynamics of a passive scalar field 8(r, t) 
advected by a velocity field v(r,t) is described by the 
following partial differential equation (written in space 
dimension d) : 

dtO(r, t) + v(r, t) ■ V B{r, t) = nV 2 6{r, t) + f(r, t), (1) 

where f(v,t) is an external source (forcing) of scalar and 
k is the molecular diffusivity. In all that follows we shall 
assume that 0(r, t) = at some distant time in the past 
t = —T (eventually, we shall let T — > oo). 

In this section the advecting velocity v and the forcing 
/ can be either deterministic or random. In the latter 
case, no particular assumption is made regarding their 
statistical properties. 

In order to illustrate the basic idea of the Lagrangian 
strategy, let us first set n — 0. We may then integrate 
(|l|) along its characteristics, the Lagrangian trajectories 
of tracer particles, to obtain 



*) = / f(a(s;r,t),s)ds, 



(2) 



where a(s; r, t) is the position at time s < t of the fluid 
particle which will be at position r at time t. (Two- 
time Lagrangian positions of this type were also used in 
Kraichnan's Lagrangian History Direct Interaction the- 
ory |1~4| .) This Lagrangian position, which will hence- 
forth be denoted just a(s), satisfies the ordinary differ- 
ential equation 



da(s) 
ds 



v(a(s),s), o(t) 



(3) 



For k > we shall now give a stochastic generalization 
of this Lagrangian representation. Roughly, 9 will be the 
average of a random field <f> which satisfies an advection- 
forcing equation with no diffusion term, in which the ad- 
vecting velocity is the sum of v and a suitable white-noise 
velocity generating Brownian diffusion. 
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To be more specific we need to introduce some no- 
tation. We shall use a set of n d-dimensional time- 
dependent random vectors 

wt(s) = {w i>a ; i = 1, . . . , n; a = 1, . . . , d}, (4) 

which are Gaussian, identically distributed, independent 
of each other and independent of both v and /. The time 
dependence is assumed to be white noise : 



(Wi^ a (s)Wj^(s')) w = 8rj5 a p8(s - s'). 



(5) 



The notation (-) w stands for "average over the tb,'s for a 
fixed realization of v and /". Similarly, stands for 
"average over v and / for a fixed realization of the Wi 's" . 
Unconditional averages are denoted just (•). Clearly, 



<•> = «■>»>„/ = «■>„/>, 



(6) 



The white noise, which is a random distribution, is here 
denoted by w(s) since it is the time derivative of the 
Brownian motion (or Wiener-Levy) process. In the nu- 
merical implementation we shall work with increments of 
w(s). 

We can now state the main result which is at the basis 
of our Lagrangian method. 

Let cf>i(r,t) (i — 1, . . . ,n) be the solutions of the fol- 
lowing advection-forcing equations : 

d t <t>i{r,t)+ (t>(r,t) + V&i(t)) • V&(r,t) = f(r,t). 
&(r,-T)=0, i = l,2,...,n. (7) 

For any T\, r 2 , r 3 , . . . , we have 

9(r 1 ) = {cl>(r 1 )) w , V* (8) 

0(ri)0(r a ) = (Mri)^(r 2 )) w , Vi ^ j, (9) 
0(ri)0(r a )0(r 3 ) = (&(ri)^(r 2 )<fe(r 3 )) w , 

Mi,j,k, i^j^k^i, (10) 



where all the fields 9 and <f> are evaluated at the same 
time t. 

The proof of (||) is obtained by taking the mean value 
of (Q) (only over w) and noting that 



2nWi ■ V</>, 



(11) 



(This is a standard result for linear stochastic equations 
having white-noise coefficients; it is derived using Gaus- 
sian integration by parts (see Ref. jl3|, Chap. 4). For 
similar derivations see Refs. §,111,11110 

To prove (ji|), we first derive from (m), assuming i ^ j, 

d t (&(ri)<k(r 2 )) + [«(n) • Vi + v(r 2 ) • V 2 
+V2K,(wi ■ Vi + wj ■ V2)]^»(ri)^j-(r 2 ) = 
/(ri)^(r a )+/(r 2 )^(ri), (12) 



where Vi and V 2 stand for Vr 1 and Vr 2 - We then 
average ([l2]) (over w) and use a relation similar to ( [ll"|) 



(13) 



^(ibi ■ Vi + w 2 ■ V 2 ) (t)i{r 1 )4> J {r 2 ) 
- K N\ + V%) (&( ri )^(r 2 )) w . 



(Notice that cross terms involving Vi • V 2 disappear be- 
cause of the independence of Wi and Wj.) The higher 
order equations in ( |l0|) and following are proved simi- 
larly. 

Because of the absence of a diffusion operator in (Q), 
its solution has an obvious Lagrangian representation 



<t>i{r,t) 



f(a>i(s),s) ds, 



where 



ddi(s) 



ds 



= v (oi(s), s) + V2nWi(s), 



Oi(t) = r. 



(14) 



(15) 
(16) 



So far we have not used the random or deterministic 
character of v and /. In the random case, taking the 
average of (|8|)~(|l0|) over v and /, we obtain 

(«(n)) = (^i)>, V* (17) 

(0(ri)0(r 2 )) = <&(n)&(r 2 )>, Vi^j (18) 
(0(ri)0(r 2 )0(r s )) = (&(n)^(r 2 )0 fe (r 3 )) , 

Vi,j,k, i^j^k^i, (19) 



Eqs. Q, (0), @, together with @, @ and 
higher orders, constitute our Lagrangian representation 
for multiple-point moments of the passive scalar. 

Let us stress that, for moments beyond the first order, 
it is essential to use more than one white noise. Indeed, 
we could have made use of just (0) and written 



(ri)0(r 2 



i{ri))w{<i>i{r 2 )) Wl 



(20) 



(Including the case i — j.) It is however not possible 
to (v, /)-average (^) because it involves a product of 
u;-averages rather than just one average as in (^). 

In the more restricted context of the Kraichnan model, 
a functional integral representation of nth order moments 
involvin g, a s here, n white noises has already been given 
in Refs. 



III. THE LIMIT OF VANISHING MOLECULAR 
DIFFUSIVITY 

Interesting pathologies occur when we bring two or 
more Eulerian space arguments of the moments to co- 
incidence and simultaneously let k — > 0. 

This is already seen on the second-order moment 
(# 2 (ri)). From (Q) and @, we have 
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! (n)> = 



(f(a i (s),s)f(a j (s'),s')) dsds', 



(21) 



for any i =/= j. The differential equations for a^(s) and 
Oj(s) involve different white noises. If, in the limit k — > 0, 
we simply ignore the ^/2nWi terms in ([15), we find that 
all the cii(s)'s satisfy the same equation (p) and the same 
boundary condition a,i(t) — r±. It is then tempting to 
conclude that all the eij(s)'s are identical, namely are 
just the Lagrangian trajectories a(s) of the unperturbed 
v-fiow. As a consequence, we can then rewrite (|2l]) as 



(22) 



»(n))= / / (f(a(s),s)f(a(s') )S ')) dsds', 

J — T J —T 

2\ 

f(a(s), s) ds 




For a large class of random forcings / of zero mean 
the r.h.s. of @ will diverge cx T as T -> oo. This is 
for example the case when / is homogeneous, station- 
ary and short- (or delta-) correlated in time as in the 
Kraichnan model. The reason of this divergence is that, 
although the forcing has zero mean, its integral (along 
the Lagrangian trajectory) over times long compared to 
the correlation time behaves like Brownian motion (in 
the T-variable) and, thus, has a variance cx T. From 
this naive procedure we would thus conclude that, when 
T = oo, the scalar variance becomes infinite as n — ► 0. 

This conclusion is actually correct when the v-Row is 
smooth (differentiable in the space variable) : this is the 
so-called Batchelor limit which has been frequently in- 
vestigated j^,|l6|,|l8 21 1 . The conclusion however becomes 



incorrect when the u-flow is only Holder continuous, i.e. 
its spatial increments over a small distance £ vary as a 
fractional power of £ (e.g. £ x ^ in Kolmogorov 1941 tur- 
bulence). As pointed out in Ref. jl6| (see also Ref. Q), 
when v is not smooth the solution to (||) lacks unique- 
ness, so that two Lagrangian particles which end up at 
the same point r± at time t may have different past histo- 
ries. This is exemplified with the one-dimensional model 



da- 
ds 



,1/3 



-T < s <t, x(t) = e > 0, (23) 



where x is a deputy for — a,j, the separation between 
two nearby Lagrangian particles. For e > 0, the solution 
of (||) is 



x(s) 



3/2 



(24) 



If we now set e = in ( pi[ ) or consider times s such that 
It — s\ > e 2 / 3 , we obtain 



x(s) 



3/2 



This is indeed a solution of (g3|) with e = 0, but there 
is another trivial one, namely x(s) = 0. Related to this 
non-uniqueness is the fact that, when e is small the solu- 
tion given by ( p4|) becomes independent of e, namely is 
given by the non-trivial solution ( p5| ) for e = 0. 

Whenever the flow v is just Holder continuous, the 
separation of nearby Lagrangian particles proceeds in a 
similar way, becoming rapidly independent of the initial 
separation. Such a law of separation is much more ex- 
plosive than would have been obtained for a smooth flow 
with sensitive dependence of the Lagrangian trajectories. 
In the latter case we would have x(s) = e e A (*~ s ) (A > 0), 
which grows exponentially with t — s but still tends to 
zero with e. 

We propose to call this explosive growth a Richard- 
son walk, after Lewis Fry Richardson who was the first 
to experimentally observe this rapid separation in turbu- 
lent flow and who was also much interested in the role of 
non-differentiability in turbulence [p2| . It is this explo- 
sive separation which prevents the divergence of (9 2 (ri)} 
when k — > (and more generally of moments with sev- 
eral coinciding points). Indeed, as the time s moves back 
from s = t, even an infinitesimal amount of molecular 
diffusion will slightly separate, say, by an amount e, the 
Lagrangian particles di(s) and CLj(s) which coincide at 
s = t. Then, the Richardson walk will quickly bring the 
separations to values independent of e and, thus of k. 
Hence, the double integral in (^l|), which involves points 
ai(s) and 02(5) with uncorrelated forces when \s — s'\ 
is sufficiently large, may converge for T — > 00. (For it 
to actually converge more specific assumptions must be 
made about the space and time correlations of v and /, 
which are satisfied, e.g., in the Kraichnan model.) 

An alternative to introducing a small diffusivity is to 
work at k — 0, with "point splitting" . For this one re- 
places (0 2 (ri)} by (9(ri)6(r'i)), where r\ and r'\ are 
separated by a distance e. Eventually, e — > 0. In practi- 
cal numerical implementations we found that point split- 
ting works well for second-order moments but is far less 
efficient than using a small diffusivity for higher-order 
moments. 



IV. THE KRAICHNAN MODEL 



(25) 



The Kraichnan model ]6|,|1|] is an instance of the pas- 
sive scalar equation ([l]) in which the velocity and the 
forcing are Gaussian white noises in their time depen- 
dence. This ensures that the solution is a Markov 
process in the time variable and that closed moment 
equations, sometimes called "Hopf equations", can be 
written for single-time multiple-space moments such as 
{6{r\, t) . . . 8(r n , t)) . The equation for second-order mo- 
ments was published for the first time by Kraichnan |J 
and, for higher-order moments, by Shraiman and Siggia 
|p0f . Note that Hopf's work [g3| dealt with the charac- 
teristic functional of random flow; it had no white-noise 
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process and no closed equations, making the use of his 
name not so appropriate in the context of white-noise 
linear stochastic equations. The fact that closed moment 
equations exist for such problems has been known for a 
long time (see, e.g., Ref. jl4| and references therein). We 
shall not need the moment equations and shall not write 
them here (for an elementary derivation, see Ref. Jjl|). 

The precise formulation of the Kraichnan model as 
used here is the following. The velocity field v — {v a , a = 
1, . . . , d} appearing in ([j]) is incompressible, isotropic, 
Gaussian, white-noise in time; it has homogeneous incre- 
ments with power-law spatial correlations and a scaling 
exponent £ in the range < £ < 2 : 



([v a (r,t)-v a {O,0)][v p (r,t) 
26{t)D aP {r), 



«/?(0,0)]) 



where, 



A»/j(r)=r* (£ + d-l)5 a(i -£ 



r a rp 



(26) 



(27) 



Note that since no infrared cutoff is assumed on the ve- 
locity its integral scale is infinite; this is not a problem 
since only velocity increments matter for the dynamics of 
the passive scalar. Note also that when a white-in-time 
velocity is used in (JT^) , the well known Ito-Stratonovich 
ambiguity could appear |p4| . This ambiguity is however 
absent in our particular case, on account of incompress- 
ibility. 

The random forcing / is independent of d, of zero 
mean, isotropic, Gaussian, white-noise in time and ho- 
mogeneous. Its covariance is given by : 



(f(r,t) f (0,0)) =F(r/L) S(t), 



(28) 



with F(0) > and F (r/L) decreasing rapidly for r ^> L, 
where L is the (forcing) integral scale. 

In principle, to be a correlation, the function F (r/L) 
should be of positive type, i.e. have a non-negative d- 
dimensional Fourier transform. In our numerical work 
we find it convenient to work with the step function 
©lM which is equal to unity for < r < L and to 
zero otherwise. Hence, the injection rate of passive scalar 
variance is e = F(0)/2 = 1/2. The fact that the func- 
tion is not of positive type is no problem. Indeed, let 
its Fourier transform be written E(k) — E\(k) — E2(k), 
where E(k) = E\(k) whenever E(k) > 0. Using the step 
function amounts to replacing in (P the real forcing / by 
the complex forcing fi + ify where /i and f 2 are indepen- 
dent Gaussian random functions, white-noise in time and 
chosen such that their energy spectra (in the space vari- 
able) are respectively E\(k) and E%(k). Since the passive 
scalar equation is linear, the solution may itself be writ- 
ten as 9\ + i02 where 9\ and 82 are the (independent) 
solutions of the passive scalar equations with respective 
forcing terms f± and f 2 . Using the universality with re- 
spect to the functional form of the forcing |25[ , it is then 
easily shown that the scaling laws for the passive scalar 
structure functions are the same as for real forcing. 



We shall be interested in the passive scalar structure 
functions of even order 2n (odd order ones vanish by 
symmetry), defined as 



S 2n (r;L) = ((9(r)-9(0)f 



(29) 



From Ref. g (see also Ref. Q) we know that, for L ^> 
r 3> r\ ~ the second-order structure function is 

given by 



S 2 (r ] L) = C 2 er^\ ( 2 = 2 - £, 



(30) 



where C2 is a dimensionless numerical constant. If there 
were no anomalies, we would have, for n > 1, 



S2n(r;L) = C 2 ne n r n ^ 



(31) 



Note that (|3(]) and (|T]) do not involve the integral scale 
L. Actually, for n > 1, we have anomalous scaling with 
S2n(r; L) oc r^ 2n and C, 2n < nCp- More precisely, we have 



S 2 n(r;L) =C 2n E n r 



/■anom „ /■ /• 

Qn = n Q2 ~ C2n, 



L 



cir n 



(32) 
(33) 



where the structure function now displays a dependence 
on the integral scale L. Our strategy will be to measure 
the dependence of S2n{i", L) on L while the separation r 
and the injection rate e are kept fixed and, thereby, to 
have a direct measurement of the anomaly Clri° m - 

Let us show that, in principle, this can be done 
by the Lagrangian method of Section [n] using 2n 
tracer (Lagrangian) particles whose trajectories satisfy 
(|l5|). The structure function of order 2n can be writ- 
ten as a linear combination of moments of the form 
(9(ri)9(r2) ■ ■ ■ 9(r 2n )), where p of the points are at loca- 
tion r and 2n — p at location (p — 0, . . . , 2n). Because 
of the symmetries of the problem, p and 2n — p give the 
same contribution. Thus, we need to work only with the 
n + 1 configurations corresponding to p = 0, . . . , n. For 
example, we have 

S 4 (r; L) = 2 (9 4 (r)) - 8 (9 3 (r)9(0)) + 6 (9 2 {r)9 2 {Q)) . 

(34) 

Let us first consider the case of the two-point function 
(second-order moment). Using ((TJ), ( |l8| ) and the inde- 
pendence of v and /, we have 



<0(ri)0(r 2 )) - 
t pt 

J (f{ai(si),si)f(a 2 (s 2 ),S2)) f dsids2^ ■ (35) 



Here, (•) / is an average over the forcing and (-) vw denotes 
averaging over the velocity and the ti;'s, and a\(si) and 
12(^2) satisfy (jl5|) with the "final" conditions cti(i) = 
riand d2(t) — r2, respectively. 
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In ( |35| ) the averaging over / can be carried out explic- 
itly using (|2^). With our step-function choice for F, we 
obtain 

{6(r 1 )6{r 2 )) = {T 12 {L)) v , (36) 

where 

Tia(£)= / 0i(|oi(a) - a 3 (s)|) ds (37) 

is the amount of time that two tracer particles arriving 
at Ti and r 2 and moving backwards in time spent with 
their mutual distance |ai(s) — 02(a) | < L. Whether the 
particles move backwards or forward in time is actually 
irrelevant for the Kraichnan model since the velocity, be- 
ing Gaussian, is invariant under reversal. 

For the four-point function, we proceed similarly and 
use the Wick rules to write fourth-order moments of / 
in terms of sums of products of second-order moments, 
obtaining 

(0(ri)0(r 2 )0(r 3 )0(r 4 )) = (T 12 {L)T 3i (L)) vw 
+ (T 13 {L)T 2i (L)) vw 
+ (T U {L)T 23 (L)) VW . (38) 

Expressions similar to Eqs. (|3^) and ( pq ) are easily de- 
rived for higher order correlations. 

We see that the evaluation of structure functions and 
moments has been reduced to studying certain statisti- 
cal properties of the random time that pairs of particles 
spend with their mutual distances less than the integral 
scale L. Generally, the distance between pairs of particles 
tends to increase with the time elapsed but, occasionally, 
particles may come very close and stay so; this will be 
the source of the anomalies in the scaling. 

V. NUMERICAL IMPLEMENTATION OF THE 
LAGRANGIAN METHOD 

In Section |v| we have shown that structure functions 
of the passive scalar are expressible in terms of u-averages 
of products of factors T,j(L). For the structure function 
of order 2n, these products involve configurations of 2n 
particles, p of which end at locations r at time s = t and 
the remaining In — p at location 0. In the Kraichnan 
model v and / are stationary, so that after relaxation 
of transients, 9 also becomes stationary. We may thus 
calculate our structure functions at t = 0. Time-reversal 
invariance of the v field and of the Lagrangian equations 
allows us to run the s-time forward rather than backward. 
Also, Tjj(L) is sensitive only to differences in particle sep- 
arations, whose evolution depends only on the difference 
of the velocities at and aj. Furthermore, the v field 
has homogeneous increments. All this allows us to work 
with 2n — 1 particle separations, namely 



a.j(s) = aj(s) - 02n(s), i = 1, ...,2n— 1, (39) 
&i(t) = Ti = n- r 2n . (40) 

Using ( |l5| ) we find that the quantities Oj(s) satisfy 2n — 1 
(vector-valued) differential equations which involve the 
differences of velocities v(a.i(s),s) — v(a 2n (s),s). The 
statistical properties of the solutions remain unchanged 
if we subtract a 2n (s) from all the space arguments. We 
thus obtain the following Lagrangian equations of motion 
for the 2n — 1 particle separations : 

^l=v i (s) + V2Hw i (s) (41) 
as 

v i (s)=v(a i (s),s)-v(0,s), (42) 

Wi(s) = Wi - w 2n . (43) 

For numerical purposes (f4l| ) is discretized in time using 
the standard Euler-Ito scheme of order one half p3] 

Oi(a + As) - Oj(s) = VAs (V t + V~2k, Wi), (44) 

where As is the time step and the V^'s and the WVs (i = 
1, . . . , 2n— 1) are (i-dimensional Gaussian random vectors 
chosen independently of each other and independently at 
each time step and having the appropriate correlations, 
which are calculated from (^|) and (pq), namely 

(Vi,c V j>0 ) = 

D a p(a,i) + D a p(aj) - D a/3 (ai - &j) , (45) 
(yV i , a W j ,p) = (l + 5i j )S a p, (46) 

where D a p{r) is defined in j27|). 

To actually generate these Gaussian random variables, 
we use the symmetry and positive definite character of 
covariance matrices like (^) and (p6|). Indeed, any such 
matrix can be factorized as a product (taking V as an 
example) p6[ 

(V ® V) = CC T , (47) 

where C is a nonsingular lower triangular matrix and C T 
is its transpose. C can be computed explicitly using the 
Cholesky decomposition method ^6|, an efficient algo- 
rithm to compute the triangular factors of positive def- 
inite matrices. It nevertheless takes 0([(2n— l)d] 3 /3) 
flops to get C from (V ® V) and this is the most time- 
consuming operation at each time step. Once C is ob- 
tained, a suitable set of variables V can be obtained by 
applying the linear transformation 

V = CM (48) 

to a set of independent unit-variance Gaussian random 
variables Afi. a coming from a standard Gaussian random 
number generator, that is with (Afi )a J\fj,f}) = SijS a p. The 
resulting variables (Vj) a then have the required covari- 
ances. 

From the cij(s)'s we obtain the quantities Ty(i) for 
all the desired values of the integral scale L, typically, a 
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geometric progression up to the maximum value L max . 
The easiest way to evaluate S2n{f",L) is to evolve si- 
multaneously the n+l configurations corresponding to 
p = 0, . . . , 7i, stopping the current realization when all 
inter-particle distances in all configurations are larger 
than some appropriate large-scale threshold L t h, to which 
we shall come back. The various moments appearing in 
the structure functions are then calculated using expres- 
sion such as ( |38| ) in which the (v, w)-averaging is done 
by the Monte-Carlo method, that is over a suitably large 
number of realizations. We note that the expressions 
for the structure functions of order higher than two in- 
volve heavy cancellations between the terms correspond- 
ing to different configurations of particles. For instance, 
the three terms appearing in the expression (34) for the 
fourth-order function, all have dominant contributions 
scaling as L 2 ( 2_ f for large L and a first subdominant 
correction scaling as L 2- ^. The true non-trivial scaling 
oc L^ 4 emerges only after cancellation of the domi- 
nant and first subdominant contributions. For small £ 
the dominant contributions are particularly large. In 
the presence of such cancellations, Monte-Carlo averag- 
ing is rather difficult since the relative errors on indi- 
vidual terms decrease only as the inverse square root of 
the number of realizations. In practice, the number of 
realizations is increased until clean non-spurious scaling 
emerges. In three dimensions, for Si(r;L), between one 
and several millions realizations (depending on the value 
of £) are required. In two dimensions even more real- 
izations are needed. For example, to achieve comparable 
quality of scaling for £ = 0.75, in three dimensions 4 x 10 6 
realizations are needed but 14 x 10 6 are needed in two 
dimensions. 

Now, some comments on the choice of parameters. 

The threshold Lth must be taken sufficiently large com- 
pared to largest integral scale of interest L max to ensure 
that the probability of returning within L max from L t h 
is negligible. But choosing an excessively large Lth is too 
demanding in computer resources. In practice, the choice 
of Lth depends both on the space dimension and on how 
far one is from the limit £ = 0. In three dimensions, it is 
enough to take L t h = 10 L max . In two dimensions there 
is a new difficulty when £ is small. At £ = the motion 
of Lagrangian particles and also of separations of pairs 
of particles is exactly two-dimensional Brownian motion. 
As it is well-known, in two dimensions, Brownian motion 
is recurrent (see, e.g., Ref. H?]]). Hence, with probability 
one a pair of particles will eventually achieve arbitrarily 
small separations. As a consequence, at £ = in two di- 
mensions, the mean square value of is infinite. For very 
small positive £ this mean square saturates but most of 
the contribution comes from scales much larger than the 
integral scale. This forces to choose extremely large val- 
ues of L t h when £ is small. In practice, for 0.6 < £ < 0.9 
we take Lth = 4 x 10 3 L max and beyond £ = 0.9 we 
take L t h = 10L max . (The range £ < 0.6 has not yet 
been explored.) In view of the accuracy of our results, 
we have verified that the use of larger values for Lth does 



not affect in any significant way the values of the scaling 
exponents. 

The molecular diffusivity k is chosen in such a way 
that r is truly in the inertial range, namely, we demand 
(i) that the dissipation scale rj ~ n 1 ^ should be much 
smaller than the separation r and (ii) that the time a 
pair of particle spends with a separation comparable to 
77, which is ~ rj 1 / k ~ 7/ 2 ~^ should be much smaller than 
the time needed for this separation to grow from r to L, 
which is ~ L 2 ~£ — r 2 ~^. The latter condition becomes 
very stringent when £ is close to 2. 

Finally, the time step As is chosen small compared to 
the diffusion time rj 2 / k at scale 77 . 



VI. RESULTS 

We now present results for structure functions up to 
fourth order. The three-dimensional results have already 
been published in Ref. ]l0| . The two-dimensional results 
are new. Some results for structure functions of order six 
have been published in Ref. p3] and shall not be repeated 
here (more advanced simulations are in progress). 
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FIG. 1. 3-D second-order structure function S2 vs L 
for £ = 0.6. Separation r = 2.7 x 10~ 2 , diffusivity 



A severe test for the Lagrangian method is provided 
by the second-order structure function S^(r;L), whose 
expression is known analytically Q[. Its behavior being 
non- anomalous, a flat scaling in L should be observed. 
The L-dependence of 5*2, measured by the Lagrangian 
method, is shown in Fig. |] for £ = 0.6 and d = 3 (all 
structure functions are plotted in log-log coordinates). 
The measured slope is 10 -3 and the error on the constant 
is 3%. (These figures are typical also for other values of £ 
studied.) We observe that, for separation r much larger 
than the integral scale L, correlations between 9{r) and 
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6(0) are very small; hence, the scaling for the second- 
order structure function is essentially given by the in- 
dependence of (6 2 ), namely L 2 ~^; the transition to the 
constant-in-L behavior around r = L is very sharp, on 
account of the step function chosen for F. 
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FIG. 2. 3-D fourth-order structure function S4 vs L for 
£ = 0.2. Separation r = 2.7 x 1CT 2 , diffusivity k = 0.247, 
number of realizations 15 x 10 6 . 
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FIG. 4. Same as in Fig. g for £ = 1.75. Parameters: 
r = 2.7 x 10~ 2 k = 10~ 9 , number of realizations 1.5 x 10 6 . 
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FIG. 5. 2-D fourth-order structure function Sa vs L for 
£ = 0.6. Parameters : r = 2.7 x 10~ 2 , K = 1.1 X 10~ 2 , number 
of realizations 5 x 10 6 . 
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FIG. 3. Same as in Fig. g for £ = 0.9. Parameters : 
r = 2.7 x 10~ 2 , K = 4.4 x 10~ 4 , number of realizations 8 x 10 6 . 
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Integral scale L 

FIG. 6. Same as in Fig. ^ for £ = 0.9. Parameters are 
r = 2.7 x 10~ 2 , k = 4.4 x 10~ 4 . To illustrate convergence, 
various numbers of realizations are shown: (1) 150 x 10 3 , (2) 
1.5 x 10 6 , (3) 3.4 x 10 6 , (4) from 4.8 x 10 6 to 7 x 10 6 (several 
curves superposed). 



contributions to ([34j), as explained in Section [v|. 

The two-dimensional case, which is numerically more 
difficult for reasons explained near the end of Section [v|, 
is shown in Figs. 0, |f and for £ = 0.6, 0.9 and 1.75, 
respectively. Fig. rcTalso shows the data obtained for var- 
ious values of the number of realizations. Note that if 
only 150 x 10 3 realizations are used, the anomaly (that is 
the slope obtained, e.g., by a least square fit) is grossly 
overestimated. 
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FIG. 7. Same as in Fig. for £ = 1.75. Parameters : 



2.7 x 10" 



10" 



number of realizations 2.4 x 10 



Figs. ||, || and || show the L-dependence of the fourth- 
order structure function in three dimensions for £ = 0.2, 
0.9 and 1.75, respectively. In each case the scaling re- 
gion (which is basically L > r) is indicated by a dashed 
straight line whose slope is the anomaly. Note that, to 
obtain a similar high-quality scaling as shown on these 
figures, a much larger number of realizations is needed for 
small £ ; this is required to permit cancellation of leading 



FIG. 8. The anomaly 2^2 — ^4 for the fourth-order structure 
function in two dimensions (stars, upper graph) and three 
dimensions (circles, lower graph). Error bars in 2-D shown 
only for £ < 1.1. The dashed line is the three-dimensional 
linear ansatz prediction d49|). 

Fig. H shows a plot of the anomaly Q nom vs £ in both 
two and three dimensions. The error bars (shown in 2- 
D only for £ < 1.1 to avoid crowding) are obtained by 
analyzing the fluctuations of local scaling exponents over 
octave ratios of values for L, a method which tends to 
overestimate errors. 

Let us now comment on the results. In three dimen- 
sions, the error bars near £ = 2 are exceedingly small 
and the data have a good fit (shown as dashed line) of 
the form Q nom — aj + 67 s / 2 with 7 = 2 — £ (the parame- 
ters are a = 0.06 and b = 1.13). This is compatible with 
an expansion in powers of ^Fj J2{| in which a term cx ^7 
is ruled out by the Holder inequality (4 < 2^2 = 27. 

Of particular significance is that, when £ is decreased 
from 2 to 0, the anomaly grows at first, achieves a max- 
imum and finally decreases. In three dimensions, where 
small-£ simulations are easier than in two dimensions, we 
have good evidence that the anomaly vanishes for £ — > 
as predicted by the perturbation theory of Gawgdzki and 
Kupiainen Q , whose leading order £| nom — 4£/5 is shown 
as a dot-dashed straight line on Fig . || Note that the 
next-order correction oc £ 2 is known pO] but the conver- 
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gence properties of the ^-series are not clear. The "linear 
ansatz" prediction for the anomaly given in Refs. [jl],f 
that is 



/■ai 
C2 



3C2 



1 



(49) 
(50) 



is consistent with our results only near £ = 1, the point 
farthest from the two limits £ = and £ = 2, which both 
have strongly nonlocal dynamics. This suggests a pos- 
sible relation between deviations from the linear ansatz 
and locality of the interactions plj . Whether this con- 
sistency for £ = 1 persists for moments of order higher 
than four is an open problem. 

The fact that the anomalies are stronger in two than 
in three dimensions is consistent with their vanishing as 
d — > oo H . The fact that the maximum anomaly occurs 
for a value of £ smaller in two than in three dimensions 
can be tentatively interpreted as follows. Near £ = 
the dynamics is dominated by the nearly ultraviolet- 
divergent eddy diffusion, whereas near £ = 2 it is dom- 
inated by the nearly infrared-divergent stretching. The 
former increases with d, but not the latter. The maxi- 
mum is achieved when these two effects balance. 



VII. CONCLUDING REMARKS 

Compared to Eulerian simulations of the partial dif- 
ferential equation ([!]) our Lagrangian method has var- 
ious advantages. When calculating moments of order 
In we do not require the complete velocity field at each 
time step but only 2n — 1 random vectors, basically, the 
set of velocity differences between the locations of La- 
grangian tracer particles which are advected by the flow 
and subject to independent Brownian diffusion. As a con- 
sequence the complexity of the computation (measured in 
number of floating point operations) grows polynomially 
rather than exponentially with the dimension d. Further- 
more, working with tracers naturally allows to measure 
the scaling of the structure functions 52 n (r; L) vs the in- 
tegral scale L of the forcing. Physically, this means that 
the injection rate of passive scalar variance (which equals 
its dissipation rate) and the separation r are kept fixed 
while the integral scale L is varied. Anomalies, that is 
discrepancies from the scaling exponents which would be 
predicted by naive dimensional analysis, are measured 
here directly through the scaling dependence on L of the 
structure functions. 

Actually, direct Eulerian simulations and the La- 
grangian method are complementary. Very high- 
resolution simulations of the sort found in Ref. (5) are 
really not practical in more than two dimensions but 
do give access to the entire Eulerian passive scalar field. 
Hence, they can and have been used to address questions 
about the geometry of the scalar field and about proba- 
bility distributions. 



Although we have presented here the numerical imple- 
mentation of the Lagrangian method only for the Kraich- 
nan model, it is clear that the general strategy presented 
in Sections || and III is applicable to a wide class of ran- 
dom flows, for example, with a finite correlation time. 

An important aspect of the results obtained for the 
Kraichnan model is that they agree with perturbation 
theory Q for £ — > 0. As is now clear from theory and 
simulations, anomalies in the Kraichnan model and other 
passive scalar problems arise from zero modes in the oper- 
ators governing the Eulerian dynamics of n-point correla- 
tion functions [ppj]32]] . It is likely that some form of zero 
modes is also responsible for anomalous scaling in nonlin- 
ear turbulence problems. An instance are the "fluxless 
solutions" to Markovian closures based on the Navier- 
Stokes equations in dimension d close to two, which have 
a power-law spectrum with an exponent depending con- 
tinuously on d p3]. Attempts to capture intermittency 
effects in three-dimensional turbulence are being made 
along similar lines (see, e.g., Ref. (34j|). The Kraichnan 
model for passive scalar intermittency puts us on a trail 
to understanding anomalous scaling in turbulence. 
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FIGURE CAPTIONS 

Fig. 1. 3-D second-order structure function S2 vs L for 
£ = 0.6. Separation r = 2.7 x 10~ 2 , diffusivity k = 
1.115 x 10~ 2 , number of realizations 4.5 x 10 6 . 

Fig. 2. 3-D fourth-order structure function S4 vs L for 
£ = 0.2. Separation r = 2.7 x 10" 2 , diffusivity k = 0.247, 
number of realizations 15 x 10 6 . 

Fig. 3 Same as in Fig. 2 for £ = 0.9. Parameters: r = 
2.7 x 10~ 2 , k = 4.4 x 10 -4 , number of realizations 8 x 10 6 . 

Fig. 4. Same as in Fig. 2 for £ = 1.75. Parameters: 
r = 2.7xl0~ 2 , k = 10~ 9 , number of realizations 1.5 x 10 6 . 
Fig. 5. 2-D fourth-order structure function S4 vs L for 
£ = 0.6. Parameters: r = 2.7 x 10~ 2 , k = 1.1 x 10~ 2 , 
number of realizations 5 x 10 6 . 

Fig. 6. Same as in Fig. 5 for £ = 0.9. Parameters are 
r = 2.7 x 10~ 2 , k = 4.4 x 10~ 4 . To illustrate convergence, 
various numbers of realizations are shown : (1) 150 x 10 3 , 
(2) 1.5 x 10 6 , (3) 3.4 x 10 6 , (4) from 4.8 x 10 6 to 7 x 10 6 
(several curves superposed). 

Fig. 7. Same as in Fig. 5 for £ = 1.75. Parameters: 
r = 2.7x 10~ 2 , k = 10~ 9 , number of realizations 2.4x 10 6 . 

Fig. 8. The anomaly 2(2 — Ct f° r the fourth-order struc- 
ture function in two dimensions (stars, upper graph) and 
three dimensions (circles, lower graph). Error bars in 
2-D shown only for £ < 1.1. The dashed line is the three- 
dimensional linear ansatz prediction (49). 
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